/*----------------------------------------------------
Class: plant_2DSOS_RDMD
DESCRIPTION:
  This class implements random deposition(RD) with migration(M)
  and desorption (D). The lattice is 2D solid-on-solid(SOS).

REFERENCE/BASED ON:
  n_sos_2D_L200_T650.c by Gangshi and Jianqiao
  sos_2d_L100_T700.c

INPUT:
 * T: Temperature, K
 * W: deposition rate, layer/s
OUTPUT:
 * z1~z4: covariance of the four modes
 * R2: roughness square, layer^2
 * M2: RMS slope square, dimensionless
 * H: average height, layer
STATE:
 * lattice: a 2D array save the height profile.
OPTIONS:
 *

AUTHOR: Xinyu Zhang (zxy1256@gmail.com)
----------------------------------------------------*/
#ifndef PLANT_2DSOS_RDMD_H
#define PLANT_2DSOS_RDMD_H

#include "../Signal.h"

const double EV2KB = 11608.696;
const double KB2H  = 2.08144796e10;

class plant_2DSOS_RDMD:public block
{
private:
	// Options
	int    LatticeSize;
	int    mode;
	double E_D;
	double E_S;
	double E_N;
	int    FLAG_CAL_MODE;
	
	// Inputs
	double T;
	double W;

	// States
	int**  lattice;
	int**  lattice_stat;
	double *p_migr;
	double *p_desorp;
	double ads_ttl,migr_ttl,desorp_ttl,rate_ttl;
	double MIGR_RT[5];
	double DESORP_RT[5];
	double k_m2d;

	// Outputs
	double *cov_z1;
	double *cov_z2;
	double *cov_z3;
	double *cov_z4;
	double h;
	double r2;
	double m2;

	// Others
	double **ef_cos_table; // A precalculated table for the cos part in eigenfunctions
	double **ef_sin_table; // A precalculated table for the sin part in eigenfunctions

	void add_remove_atom(int x, int y, int num);
	void adsorption();
	void migration();
	void desorption();
	double OneMCTime(double r_d);
	void migr_site(int *migr_x,int *migr_y);
	void desorp_site(int *desorp_x,int *migr_y);
	void set_ads_rate(double W_set);
	void set_temperature(double T_set);
	void GetNeighbors(int x,int y, int *xpos,int *ypos);	
	void init_migr_rate();
	void init_desorp_rate();
	void init_p_migr();
	void init_p_desorp();
	void init_local_stat(int x, int y);
	void h2mode();
	void h2mode_fftw();

	inline int index2D_h(int i,int j){return i*LatticeSize+j;};
	inline int index2D_mode(int m,int n){return m*(mode+1)+n;};

	// Eigen-functions
	inline double phi1mn(int m,int n,int ix,int iy){
		//return 2.0*sin(2*m*pi/LatticeSize*ix)*sin(2*n*pi/LatticeSize*iy)/pi;
		return ef_sin_table[m][ix]*ef_cos_table[n][iy];
	};

	inline double phi2mn(int m,int n,int ix,int iy){
		/*double coeff;
		if((m == 0 || m*2==LatticeSize)&&(n==0 || n*2==LatticeSize)){
			coeff = 1.0/pi;
		}
		else if(m!=0 && 2*m!=LatticeSize && n!=0 && 2*n!=LatticeSize){
			coeff = 2.0/pi;
		}
		else{
			coeff = SQRT_TWO/pi;
		}
		return coeff*cos(2*m*pi/LatticeSize*ix)*cos(2*n*pi/LatticeSize*iy);*/
		return ef_cos_table[m][ix]*ef_sin_table[n][iy];
	};

	inline double phi3mn(int m,int n,int ix,int iy){
		/*double coeff;
		if(n==0||n*2 == LatticeSize){
			coeff = SQRT_TWO/pi;
		}else{
			coeff = 2.0/pi;
		}
		return coeff*sin(2*m*pi/LatticeSize*ix)*cos(2*n*pi/LatticeSize*iy);*/
		return ef_sin_table[m][ix]*ef_sin_table[n][iy];
	};

	inline double phi4mn(int m,int n,int ix,int iy){
		/*double coeff;
		if(m==0||2*m==LatticeSize){
			coeff = SQRT_TWO/pi;
		}else{
			coeff = 2.0/pi;
		}
		return coeff*cos(2*m*pi/LatticeSize*ix)*sin(2*n*pi/LatticeSize*iy);*/
		return ef_cos_table[m][ix]*ef_sin_table[n][iy];
	};

	// The following declaration is to prevent compiler default implementation.
	plant_2DSOS_RDMD(const plant_2DSOS_RDMD&);
	plant_2DSOS_RDMD& operator=(const plant_2DSOS_RDMD&);
public:
	plant_2DSOS_RDMD(double idt,int iLx=100,int iMode = 50,double iE_D = 1.8,double iE_S = 1.58,double iE_N=0.28);
	plant_2DSOS_RDMD(const string FileName);
	~plant_2DSOS_RDMD();
	void update(double sysTime,double iT,double iW);
	virtual void update(double sysTime);
	void   reset();
	void   setNoMode(){FLAG_CAL_MODE = 0;};
	void   getSnapshot(const string FileName);
	double getM2();
	double getR2();
	double getH();
	double getTime();
};

double ranO(int *idum);
double RandRand();
double RandN1();
#endif
